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ABSTRACT 

We present a model for star formation and supernova feedback that describes the 
multi-phase structure of star forming gas on scales that are typically not resolved 
in cosmological simulations. Our approach includes radiative heating and cooling, 
the growth of cold clouds embedded in an ambient hot medium, star formation in 
these clouds, feedback from supernovae in the form of thermal heating and cloud 
evaporation, galactic winds and outflows, and metal enrichment. Implemented using 
smoothed particle hydrodynamics (SPH), our scheme is a significantly modified and 
extended version of the grid-based method of Yepes et al. (1997), and enables us to 
achieve high dynamic range in simulations of structure formation. 

We discuss properties of the feedback model in detail and show that it predicts 
a self-regulated, quiescent mode of star formation, which, in particular, stabilises the 
star forming gaseous layers of disk galaxies. The parameterisation of this mode can be 
reduced to a single free quantity which determines the overall timescale for star forma- 
tion. We fix this parameter numerically to match the observed rates of star formation 
in local disk galaxies. When normalised in this manner, cosmological simulations em- 
ploying our model nevertheless overproduce the observed cosmic abundance of stellar 
material. We are thus motivated to extend our feedback model to include galactic 
winds associated with star formation. 

Using small-scale simulations of individual star-forming disk galaxies, we show 
that these winds produce either galactic fountains or outflows, depending on the depth 
of the gravitational potential. In low-mass haloes, winds can greatly suppress the over- 
all efficiency of star formation. When incorporated into cosmological simulations, our 
combined model for star formation and winds predicts a cosmic star formation density 
that is consistent with observations, provided that the winds are sufficiently energetic. 
Moreover, outflows from galaxies in these simulations drive chemical enrichment of the 
intcrgalactic medium, in principle accounting for the presence of metals in the Lyman 
alpha forest. 
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1 INTRODUCTION 

The past two decades have witnessed the emergence of a suc- 
cessful class of theoretical models to explain the formation of 
structure in the Universe. In these "hierarchical" scenarios, 
cold dark matter (CDM) forms the dominant mass compo- 
nent, and structure grows via gravitational instability from 
primordial perturbations seeded by inflation. Galaxy for- 
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mation proceeds "bottom-up" in this picture, with smaller 
galaxies forming first, and larger objects being assembled 
through merging of less massive building blocks. 

Numerical simulations have become an important tool 
for exploring the detailed predictions of these models. For 
example, the clustering properties of matter on large scales, 
where gasdynamical forces are negligible compared to grav- 
ity, can now be computed with high precision. However, 
extending these studies to understand the formation of lu- 
minous components of galaxies is considerably more diffi- 
cult. Cosmological simulations that attempt to include the 
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essential physical processes associated with galaxy forma- 
tion, namely radiative cooling and star formation, invariably 
encounter formidable difficulties. For this reason, the most 
popular approach to date for describing hierarchical galaxy 
formation has been the so-called "semi-analytic" technique 
(e.g. White & Frenk, 1991; Lacey & Silk, 1991; Kauffmann 
et al., 1993, 1994; Cole et al., 1994), in which highly sim- 
plified prescriptions for the baryonic physics are combined 
with Monte-Carlo methods to construct merging history 
trees. In recent times, these Monte-Carlo trees have been 
replaced with merging histories obtained directly from high- 
resolution N-body simulations (Kauffmann et al., 1999a,b; 
Springel et al., 2001a). 

Some impressive successes in direct numerical simu- 
lation of cosmological hydrodynamics have been obtained 
since this method was pioneered more than a decade ago 
(Evrard, 1988, 1990; Hernquist, 1989; Hernquist & Katz, 
1989; Navarro & Benz, 1991; Barnes & Hernquist, 1991; 
Katz & Gunn, 1991; Hiotelis & Voglis, 1991; Thomas & 
Couchman, 1992; Katz et al., 1992; Cen & Ostriker, 1993). 
For example, these studies improved our understanding of 
quasar absorption line systems (e.g. Cen et al., 1994; Zhang 
et al., 1995; Hernquist et al., 1996), the intergalactic medium 
(e.g. Cen & Ostriker, 1993, 2000; Katz et al., 1996, 1999; 
Weinberg et al., 1997, 2002; Steinmetz & Miiller, 1995; 
Bryan et al., 1994; Bryan & Norman, 1998; Blanton et al., 
1999; Pearce et al., 1999, 2001), and galaxy formation (e.g. 
Katz & Gunn, 1991; Navarro & White, 1994a; Steinmetz & 
Miiller, 1995; Mihos & Hernquist, 1996; Walker et al., 1996; 
Navarro & Steinmetz, 1997; Steinmetz & Navarro, 1999). 

However, the numerical modelling of key processes re- 
lated to star formation and associated feedback by super- 
novae and stellar winds has remained problematic. Because 
radiative cooling is particularly efficient in dense, low-mass 
haloes at high redshift, simulations are confronted with an 
overcooling problem, which is typically manifested through 
an overproduction of stars. Semi-analytic treatments indi- 
cate that a physical handling of feedback processes is re- 
quired to regulate star formation and reduce the fraction 
of baryons that collapse in the centres of haloes. Unfortu- 
nately, it is not clear how best to incorporate feedback pro- 
cesses into simulations of galaxy formation, although they 
appear to be crucially important to a successful description 
of hierarchical CDM universes. 

Here, we present a model for star formation and feed- 
back that goes beyond a purely phenomenological descrip- 
tion. In most attempts to simulate galaxy formation, the gas 
within galaxies is treated as a single phase medium and is 
converted into stars on some characteristic timescale, which 
is usually related to the local dynamical time. Often, auxil- 
iary criteria requiring that the flow be locally converging or 
Jeans unstable are also invoked to restrict the eligibility of 
gas fluid elements for star formation. Based on an assumed 
stellar initial mass function (IMF) , it is then straightforward 
to work out the expected rate of energy input due to associ- 
ated supernovae. However, it is not clear how this "feedback 
energy" should be deposited into the gas. Simply augment- 
ing the thermal reservoir of the gas has little effect, because 
the cooling time in the dense single-phase medium in these 
simulations is so short that the energy is promptly radiated 
away. Consequently, pure thermal feedback in this manner 
fails to solve the overcooling problem. 



Other approaches have been attempted to prevent this 
rapid loss of feedback energy, thereby providing a back- 
reaction on star formation (see Kay et al., 2002, for a re- 
cent comparison of the most common of these methods). For 
example, Thackcr & Couchman (2000) temporarily lower 
the densities of SPH particles receiving feedback energy to 
prevent immediate energy loss. Alternatively, such particles 
have been treated adiabatically for a brief period of time 
after acquiring feedback energy (Gerritsen & Icke, 1997). 
Although crude, these methods show some success in regu- 
lating star formation. For example, Thacker & Couchman 
(2001) have recently found that the angular momentum 
problem in the formation of disk galaxies can be alleviated 
through the combined effects of a smoothed deposition of 
feedback energy with a delay of cooling by 30 Myr. 

Another method for depositing the energy is through ki- 
netic feedback, where radial momentum kicks are imparted 
to particles surrounding sites of star formation (e.g. Navarro 
& White, 1993; Mihos & Hernquist, 1994), but the fraction 
of feedback energy that must be put into this channel de- 
pends on numerical details of the simulations. Also, the size 
of the affected regions depends on spatial resolution, and is 
in general much larger than the physical scale of individual 
supernova remnants. 

Sommer-Larsen et al. (1999) experimented with yet 
more extreme versions of feedback in order to investigate the 
effects of galactic outflows. For example, in their "blow-out" 
simulations, they identified all dense, cold gas at z = 2.4 
and instantly removed it and redistributed it uniformly in 
spherical shells around the centres of the gas clumps. While 
this is clearly a highly artificial model, it demonstrated that 
outflows can have a large impact on galaxy formation. In 
particular, the loss of angular momentum by the gaseous 
component was greatly reduced in their simulations, allevi- 
ating the apparent discrepancy between observed and sim- 
ulated values of the specific angular momentum of galactic 
disks. Scannapieco et al. (2001) used a similar gas-removal 
technique to model galactic winds at high redshift. 

A disadvantage common to many of these recipes is 
their explicit or implicit dependence on numerical details 
of the simulations. For example, if a criterion for star for- 
mation explicitly involves the mass of a particle, it is unclear 
in what sense numerical convergence for the star formation 
rate of a specific object can be expected if the simulation is 
repeated with a different mass resolution. Instead, the free 
parameters describing star formation may need to be recali- 
brated whenever the properties of the simulation are altered. 
Springel (2000) tried to overcome this problem by developing 
a feedback method based on an effective model for the star- 
forming interstellar medium (ISM), where turbulent motions 
on unresolved scales were postulated in an ad-hoc fashion to 
provide pressure support for the gas, thereby regulating star 
formation. In essence, feedback energy was only slowly ther- 
malised in this method, thereby delaying the dissipation of 
feedback energy in a similar way to some of the models noted 
above. 

Because the Lagrangian nature of smoothed particle hy- 
drodynamics makes it relatively easy to obtain the high spa- 
tial dynamic range needed for simulations of galaxy forma- 
tion, it has so far been the most popular numerical technique 
for studying this problem. However, star formation and feed- 
back have also been incorporated into hydrodynamical mesh 



© 0000 RAS, MNRAS 000, 000-000 



Cosmological SPH simulations: A hybrid multi-phase model for star formation 3 



codes (e.g. Cen & Ostriker, 1993; Gnedin, 1998), but their 
limited spatial resolution is a serious concern, because star 
formation occurs at substantially higher overdensities than 
can be resolved with a fixed grid in cosmological volumes. 
This is quite different when high-resolution Eulerian codes 
are applied to simulations of isolated galaxies, as in the work 
of Mac Low & Ferrara (1999), who carried out detailed stud- 
ies of starbursts in dwarf galaxies. Eulerian codes have also 
been employed to directly resolve the multi-phase structure 
of the ISM, yet state-of-the-art simulations are either re- 
stricted to two dimensions (Wada & Norman, 1999), or can 
presently only resolve regions of size ~100 pc or so in three 
dimensions (Wada, 2001). 

Yepes et al. (1997) adopted an intermediate approach 
by developing a sub-grid model of a two-phase ISM and im- 
plementing it in an Eulerian code. Their spatial resolution 
was consequently relatively poor, despite restricting them- 
selves to rather small cosmological volumes, but together 
with subsequent work (Elizondo et al., 1999a, b; Ascasibar 
et al., 2002) they were able to demonstrate a number of in- 
teresting properties of this scheme, which is physically bet- 
ter motivated than most of the models currently employed in 
cosmological SPH simulations. Because of this, we will base 
our analysis in part on their approach and on a first im- 
plementation of it in SPH by Huffman & Pharasyn (1999), 
yet we will significantly extend and modify the model. Our 
specific aim is to formulate a description of feedback which 
is physically motivated, numerically well specified, and suit- 
able for simulating galaxy formation within large cosmo- 
logical volumes. In what follows, we show that our model 
leads to self-regulation of star formation, and describe how 
it can be normalised by analysing its properties. Under the 
assumption that the normalisation obtained from observa- 
tions of disk galaxies in the present Universe holds at all 
redshifts, we find that the resulting cosmic star history leads 
to an overprediction of the luminosity density of the uni- 
verse. We then argue that galactic winds are a plausible 
mechanism for eliminating this discrepancy, and they simul- 
taneously help to account for other observational data, like 
the metals found in the low-density intergalactic medium 
(IGM) , and describe a phenomenological strategy for includ- 
ing these winds in cosmological simulations. 

This paper is organised as follows. In Section 2, we 
summarise our new hybrid multi-phase model for quiescent 
star formation in detail, and analyse some of its basic prop- 
erties. In Section 3, we discuss physical and observational 
constraints on the values of the model parameters. In Sec- 
tion 4, we extend the model by including galactic winds. We 
then describe our numerical implementation of the method 
in a parallel TreeSPH code for cosmic structure formation 
in Section 5. We present numerical results for the forma- 
tion of disk galaxies in isolation, and for small simulations 
of cosmic structure formation in Section 6. Finally, we give 
conclusions in Section 7. 



2 A MULTI-PHASE MODEL FOR QUIESCENT 
STAR FORMATION 

We term our treatment of star formation and feedback a 
"hybrid" method because it does not attempt to explicitly 
resolve the spatial multi-phase structure of the ISM on small 



scales, but rather makes the assumption that crucial aspects 
of the global dynamical behaviour of the ISM can be char- 
acterised by an effective "sub-resolution" model that uses 
only spatially averaged properties to describe the medium. 
In essence, by adopting a statistical formulation, we seek to 
account for the impact of unresolved physics on scales that 
are resolved. 

In our hybrid approach, each SPH fluid element rep- 
resents a region in the ISM, whose properties are obtained 
from a suitable coarse-graining procedure. This is analogous 
to the N-body representation of collisionless fluids, but here 
the structure of the unresolved matter is more complex. In 
particular, we picture the medium as a fluid comprised of 
condensed clouds in pressure equilibrium with an ambient 
hot gas. The clouds supply the material available for star 
formation. For these hybrid-particles, the equations of hy- 
drodynamics are only followed for the ambient gas. The cold 
clouds are subject to gravity, add inertia, and participate in 
mass and energy exchange processes with the ambient gas 
phase. These processes are computed on a particle by par- 
ticle basis in terms of simple differential equations using a 
specific model for the physics of the ISM. Here, we attempt 
to incorporate some of the key aspects of the theoretical 
picture of the ISM outlined by McKee & Ostriker (1977). 

In the following, p h denotes the local density of the hot 
ambient gas, p c is the density of cold clouds, the density 
of stars, and p = pn + p c is the total gas density. Individ- 
ual molecular clouds and stars cannot be resolved, thus p c 
and p* represent averages over small regions of the ISM. We 
implicitly assume that such a coarse-graining procedure has 
been carried out, and we will formulate the interactions be- 
tween the phases assuming that the regions used to define 
the averages are of constant volume. In our model, the aver- 
age thermal energy per unit volume of the gas can then be 
written as e = phUh + PcU c , where uh and u c are the energy 
per unit mass of the hot and cold components, respectively. 

We model three basic processes that drive mass ex- 
change between the phases. These are star formation, cloud 
evaporation due to supernovae, and cloud growth due to 
cooling. We discuss these processes in turn, focusing first on 
self-regulated, "quiescent" star formation. Later, we will ex- 
tend the model to include additional processes leading to the 
development of galactic winds. For simplicity of presenta- 
tion, we will omit adiabatic terms in the following definition 
of the model. 

We assume that star formation converts cold clouds into 
stars on a characteristic timescale t+, and that a mass frac- 
tion /3 of these stars are short-lived and instantly die as 
supernovae. This can be described by: 

In principal, there is a time delay between star formation 
events and associated supernovae equal to the approximate 
lifetime of massive stars (~3x 10 7 years). However, for the 
quiescent mode of star formation, self-regulation is estab- 
lished so quickly that this time delay can be neglected. 

Star formation therefore depletes the reservoir of cold 
clouds at the rate pc/t+, and leads to an increase in the 
mass of the ambient phase as /?p c /i*, because we assume 
that ejecta from supernovae are returned as hot gas. The 
parameter /3 is the mass fraction of massive stars (> 8M ) 
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formed for each initial population of stars and hence depends 
on the adopted stellar initial mass function (IMF). For a 
Salpeter IMF (1955) with slope —1.35 and upper and lower 
limits of 40 Mq and 0.1 Mq, respectively, it has the value 
P = 0.106. We will typically adopt /3 = 0.1, and note that 
our results are not particularly sensitive to this choice. 

In addition to returning gas (enriched with metals) to 
the ambient phase of the ISM, supernovae also release en- 
ergy. The precise amount of this energy depends on the IMF. 
For the canonical value of 10 51 ergs per supernova, we ex- 
pect an average return of esn = 4 x 10 4S ergsMg 1 for each 
solar mass in stars formed for the IMF adopted here. The 
heating rate due to supernovae is hence 

d, , dp, „ p c (2) 



dt 



(phUh) 



dp* a Pc 
= e SN-T— = PUSN — , 
SN dt t* 



where usn = (1 — P)P 1 £SN may be expressed in terms of an 
equivalent "supernova temperature" Tsn = 2/iusn/(3/c) ~ 
10 8 K. 

In our formulation, we assume that the feedback en- 
ergy from supernovae directly heats the ambient hot phase. 
In addition, we suppose that cold clouds are evaporated in- 
side the hot bubbles of exploding supernovae, essentially by 
thermal conduction, thereby returning material from con- 
densed clouds to the ambient gas. We take the total mass of 
clouds that are evaporated to be proportional to the mass 
in supernovae themselves, viz. 

dp c 



dt 



(3) 



The efficiency A of the evaporation process is expected to be 
a function of the local environment. For simplicity, we will 
only take the expected theoretical dependence on density, 
A (X p~ 4/r ', into account (McKce & Ostriker, 1977). We 
will not try to estimate the evaporation efficiency from first 
principles, but rather treat its normalisation as a parameter. 
As we argue below, the range of permissible values for A is 
limited once we require a plausible temperature structure 
for the ISM. 

Finally, we invoke a process by which cold clouds come 
into existence and grow. We assume that a thermal insta- 
bility operates in the region of coexistence between the cold 
and hot phases, leading to mass exchange between ambient 
gas and cold clouds (Field, 1965; Field et al., 1969; Begel- 
mann & McKee, 1990). In particular, we argue that radiative 
energy loss by the hot gas leads to a corresponding growth 
of the clouds, i.e. the energy radiated cools gas from the 
temperature of the hot phase to that of the cold phase and 
thus gives rise to a growth in the mass of cold clouds. This 
mass flux is thus: 



dp c 
dt 



dp h 
dt 



1 



Uh 



-A net (p h ,u h ). 



(4) 



Note that we assume that the temperatures of and total 
volumes occupied by the hot and cold phases remain con- 
stant during cloud growth. We compute the cooling function 
A not from the radiative processes appropriate for a primor- 
dial plasma of hydrogen and helium essentially as described 
by Katz et al. (1996). The abundances of the various ion- 
isation states of H and He are computed explicitly, under 
the assumption of ionisation equilibrium in the presence of 
an external UV background field, which we here take to be 
a modified Haardt & Madau (1996) spectrum with reioni- 



sation taking place at redshift z ~ 6 (for details, see Dave 
et al., 1999). 

The minimum temperature the gas can reach due to 
the radiative cooling processes we include is about 10 4 K, 
at which point the gas becomes neutral, and further cooling 
would require molecular cooling or that metals be present. 
In fact, real molecular clouds will form cores that are much 
colder than 10 4 K. However, in our present statistical for- 
mulation we neglect the internal structure of the clouds. 
Their exact thermal energy content will be unimportant in 
the current model as long as Uh 3> u c , where we identify 
u h with the hot intercloud medium resulting from super- 
nova remnants at temperatures around ~ 10 5 — 10 7 K. We 
will typically assume a temperature of T c ~ 1000 K for the 
cold clouds, but note that the results do not depend on this 
choice for T c <C 10 4 K. 

In the following, we will use a factor / to differentiate 
between situations where the thermal instability is assumed 
to be operating (/ = 0), and regions where it is ineffective 
and ordinary cooling takes place (/ = 1). In general, the gas 
can be expected to exhibit thermally unstable behaviour if 
the cooling rate is a declining function of temperature. For 
a primordial mixture of gas, this is the case roughly in the 
range T w 10 5 — 10 6 K, while for metal-enriched gas the 
minimum of the cooling curve shifts to higher temperatures, 
up to a few times 10 7 K. We thus require that the tempera- 
ture of the ambient gas is above 10 s K. For the onset of the 
thermal instability itself, we chose a simple density thresh- 
old criterion, i.e. / = for p > p th , and / = 1 otherwise. 
This is motivated by the observed threshold behaviour of 
star formation (Kennicutt, 1989). Note that star formation 
will proceed only once clouds can form; i.e. when the gas 
density exceeds p t h- 

Quantitatively, the rates at which the masses of the hot 
and cold phases evolve can be written as: 



(5) 
(6) 



dp c 


= - P f-Af5 p f 


1 


-f 


A nct (ph 


Uh), 


dt 


+ — 


— u c 


dp h 
dt 




1 

Uh 


-f 

— U c 


Anot (Ph 


Uh)- 



In both of these equations, the first term on the right hand 
side describes star formation and feedback, the second cloud 
evaporation, and the third the growth of clouds due to ra- 
diative cooling of the gas. 

We now consider the energy budget of the gas, which 
we write as 



dt 



(phUh+PcUc) = -Anct(ph,U h )+P^USN-(l-P)^Uc,{7) 



where A not is the radiative cooling (or heating) function for 
the ambient hot medium. Only the ambient gas is assumed 
to be subject to radiative processes. The second term de- 
scribes the non-gravitational energy injected by exploding 
supernovae. Finally, the third term describes the loss of en- 
ergy in the gaseous phase due to material that is locked up 
in collisionless stars. The loss term arises because we as- 
sume that the material that is converted into stars is at the 
temperature of the cold clouds. 

Within the framework of our model assumptions, we 
can split the energy equation into two separate relations for 
the thermal budget of the hot and cold components: 
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-{p cUc ) = --u c 



.^ Wc+ (WK Anet; 

U u h - u c 



± ( Ph u h ) 



oPc , , \ | \ p Pc U h - fu c 
1— ("SN + «c) + A/3 — U c 

i* i* u h - u c 



(8) 



A net . (9) 



In these equations, the first term on the right hand side 
describes the effects of star formation and supernovae, the 
second accounts for cloud evaporation, and the third term 
the impact of thermal instability. 

We assume that the cold clouds are at a fixed temper- 
ature, and so u c will be treated as constant. Consequently, 
the above equations imply that the temperature of the hot 
phase will evolve according to 

P h ~JT = P^(usN + u c -Uh) - A/3^(uh-u c ) -/A nc t.(10) 

This is one of the basic equations that is integrated in the 
simulation code, augmented of course by terms coming from 
the external work done by pressure forces and viscous effects. 
In practice, we follow the evolution of the entropy of the hot 
phase in our code, rather than the thermal energy, for nu- 
merical reasons that are discussed by Springel & Hernquist 
(2002a). 

An inspection of equation (10) reveals some interesting 
properties of our model. Suppose that the thermal instability 
is operating (/ = 0), in which case the temperature of the 
hot phase changes only due to star formation and feedback. 
In fact, as long as star formation is active, the temperature 
will evolve towards 



u h = 



+ U c . 



(11) 



Deviations from this temperature decay on a timescale 



Th 



t* P h 



(12) 



Thus, provided star formation is sufficiently rapid compared 
to adiabatic heating or cooling due to gas motions, the tem- 
perature of the hot phase will be maintained at the value 
set by equation (11), which is independent of the star for- 
mation timescale i*. Since usually A 3> 1 and usn/A ^> u c 
this temperature is in practice just given by uu £± usn/A. 
For a supernova temperature of 10 8 K and A ~ 100, the self- 
regulated temperature of the hot diffuse phase of the ISM 
will thus be maintained at ~ 10 6 K. 

A further interesting consequence of our feedback model 
is that it leads to self-regulated star formation. Owing 
to evaporation, star formation reduces the density in cold 
clouds, lowering the star formation rate. On the other hand, 
a higher density of hot gas leads to an increase in the cooling 
rate, and hence to more rapid replenishing of clouds, increas- 
ing the star formation rate. In this manner, a self-regulated 
cycle of star formation is established where, in equilibrium, 
the growth of clouds is balanced by their evaporation due to 
supernova feedback. 

This condition can be seen in more detail by considering 
equation (7). In the self-regulated regime, we expect the 
effective pressure of the medium 



-Pcff = (7 - ^){phU h + PcU c ) 

to be constant in time. This condition implies 

Pc _ Anct(Ph,Mfe) 

(3u S n - (1 - P)u c ' 



(13) 



(14) 
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Figure 1. Top panel shows the temperature structure of the 
multi-phase medium as a function of baryonic overdensity (as- 
suming f2j, = 0.04). Below a density of 5 ~ 8 X 10 5 , star formation 
does not occur and the gas is treated as a single-phase fluid, with 
a temperature that is maintained close to 10 4 K by cooling and 
UV-heating. When star formation sets in, a hot ambient medium 
develops (thin solid line), and cold clouds are present at a fidu- 
cial temperature of 10 3 K (dashed). The mass- weighted effective 
temperature (thick solid line) is continuous at the onset of star 
formation, and can be interpreted in terms of an effective pres- 
sure P e ff. In the middle panel, we plot the local polytropic index 
n c ff = dlogP c ff /dlogp of the resulting effective equation of state. 
The bottom panel shows the mass fraction in cold clouds as a 
function of overdensity. 
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Using A nct (p h ,u h ) = [p h /p] 2 A nct (p, u h ), we thus obtain an 
expression for the expected ratio of masses in the cold and 
hot phases, viz. 



Pc ph 

— = — y, 

ph p 

where we have defined 
t*A net (p, Uh) 



V = 



p [/3usn - (1 - P)u c ] ' 



(15) 



(16) 



Note that y is a function only of p in the self-regulated 
regime, provided that t* and A depend only on density. We 
can then express the mass fraction 



_ Pc 

x = — 



of cold clouds as 



1 / 1 1 

2y V y 4y 2 



(17) 



(18) 



The effective pressure of the gas will then take on the value 



-Pcff = (7 - 1)P [(1 - X ) U h + XU C ] 



(19) 



where the term in square brackets is the effective mass- 
weighted "temperature" u e s of the medium. 

In Figure 1, we show the density dependence of this 
effective temperature, and of the temperatures of the hot 
and cold components in our multi-phase model of the ISM. 
Also shown is the mass fraction of gas in cold clouds, and 
the local logarithmic slope n e g of the effective equation of 
state -P e ff(p) of the star- forming medium. It is clear that the 
functional dependence of P e ff on density is particularly im- 
portant for the dynamical stability of star-forming regions. 
For n e ff > 4/3, the effective pressure can provide enough 
vertical thickening to stabilise gaseous disks against rapid 
break-up into clumps due to dynamical instabilities. 



3 SELECTING PARAMETERS 

Following McKee & Ostriker (1977), we express the density 
dependence of the supernova evaporation parameter A as 



A{p) = A (-P- 
\Pth 



-4/5 



(20) 



where p t h and Ao are parameters of our model. To spec- 
ify the star formation timescale, i*, we make the common 
assumption that this quantity is proportional to the local 
dynamical time of the gas. This plausible choice results in a 
Schmidt-type law for the dependence of the star formation 
rate on density, as observed. We thus set 



U(p) = to 



Pth 



-1/2 



(21) 



where to is an additional parameter of the model. 

Before we study the properties of our model in detail, 
we discuss how its free parameters can be constrained. To 
this point, we have introduced two parameters that depend 
on the IMF, (3 and msn, and three parameters that deter- 
mine the regulation of the multi-phase medium due to star 
formation; these are p th , Ao, and to- 



For the purposes of this work, we neglect uncertainties 
in the IMF and treat (3 and msn as known constants. How- 
ever, the other three parameters are crucially important to 
the behaviour of the model. We proceed by first requiring 
that at the onset of thermal instability, the equilibrium tem- 
perature of the hot medium is such that the thermal insta- 
bility in fact becomes operative. This means that the cooling 
function should start to fall at this temperature, which hap- 
pens at ~ 10 5 K. We thus require that 



Tsn/Aq = 10 5 K, 



(22) 



which fixes Ao to a value of Ao ~ 1000. 

Next, we argue that the effective pressure of the gas 
should be a continuous function of density at the onset of 
the regime of self-regulated star formation. The gas just be- 
low the threshold will have cooled down to 10 4 K, where it 
becomes neutral. Further cooling could happen only due to 
less efficient molecular cooling processes that we ignore. (For 
now, we also neglect cooling due to metals, which could how- 
ever become significant once the gas becomes chemically en- 
riched from star formation.) Hence we impose the condition 
Meff(pth) = 114, where ua is the specific energy corresponding 
to a temperature of 10 4 K. Based on the equations for the 
self-regulated regime, we thus have 

/3«sn - (1 - P)u c 



Pth 



Xth 



t°A{u S N/Ao) 



(23) 



(1 - Xth) 

which sets the density threshold pth for a given a value of to ■ 
Here, a; t h = 1 + (A + 1) (u c — M4) /msn — 1 — A U4 /msn is the 
mass fraction in clouds at the threshold, and A is defined by 
A(p,u) — A nct (p,u)/p 2 , such that A looses its dependence 
on density for temperatures high compared to 10 4 K. 

We have thus reduced the description of our model to 
one free parameter, the star formation timescale to- It is 
clear that this parameter sets the overall gas consumption 
timescale, which is, in principle, well-constrained by obser- 
vations of the efficiency of star formation. 

Observationally, there appears to be a tight correlation 
between disk-averaged measurements of the star formation 
rate per unit area and the gas surface density, a "global 
Schmidt-law", given by Kennicutt (1998) as 



Sspr = (2.5 ±0.7) x 10" 



M pc- 



Mr. 



yr kpc : 



.(24) 



This correlation is valid from typical disk-averaged gas den- 
sities of 10M Q pc~ 2 in normal spirals to gas densities as 
high as 10 5 M Q pc~ 2 in the central regions of starbursting 
galaxies. Roughly the same law also appears to hold locally 
for azimuthally averaged gas densities and star formation 
rates, although perhaps with a slightly smaller amplitude 
than the one given by equation (24). In addition, there is 
a clear threshold behaviour for star formation, with Esfr 
very rapidly falling at densities below ~ 10 M Q pc~ 2 (Ken- 
nicutt, 1989, 1998; Martin & Kennicutt, 2001). It has been 
suggested that the threshold might be related to the onset 
of gravitational instabilities in the disk, but presently the 
evidence for this argument remains inconclusive. 

By modelling individual galaxies, we will directly try to 
fit the Schmidt-law of equation (24), using however a slope 
of 1.5, which prompts us to lower the normalising coefficient 
in front by a factor of 2. This correction results when chang- 
ing the best-fit slope to 1.5 and simultaneously requiring 
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Figure 2. Star formation rate per unit area versus gas surface 
density in a star-forming disk of gas. The dashed inclined line 
shows the Kennicutt law of equation (25), with the vertical line 
drawn to indicate the observational cut-off in the star formation 
rate. The solid and dotted lines have been computed numerically 
by solving the equations of hydrodynamic equilibrium together 
with our multi-phase model for self-gravitating sheets of gas. The 
solid line shows the result for = 2.1 Gyr, the lower dotted 
line is obtained for = 8.4, Gyr, and the upper dotted line for 
t* = 0.53 Gyr. 

the star formation rate to remain unchanged at the inter- 
mediate density value of 10 3 Mopc~ 2 . The resulting local 
Schmidt-law then provides a good fit to the azimuthal data 
of 21 spirals presented in Fig. 3 of Kennicutt (1998). The 
observational constraint we thus try to match can also be 
expressed as a gas consumption timescale, taking the form 

tsFR = ^ = 3.2 Gyr ( _ \ ^ . (25) 

Ssfr ^lOMopc 2 J 

At the threshold for the onset of star formation, this thus 
indicates a timescale of about 3.2 Gyr, becoming shorter 
towards higher densities. This compares well with the cited 
median gas consumption timescales of 2.1 Gyr (Kennicutt, 
1998) and 2.4 Gyr (Rownd & Young, 1999). 

In Figure 2, we show the relation between star forma- 
tion density and gas surface density predicted by our multi- 
phase model for various choices of ij. The curves have been 
computed by solving the equations of hydrostatic equilib- 
rium for self-gravitating layers of gas with varying surface 
density. Clearly, the amplitude of the star formation rate is 
very sensitive to the value of to, with to = 2.1 Gyr providing 
a good fit to the Kennicutt law. 

Once the amplitude of star formation has been matched 
by adjusting to, the slope and the cut-off obtained can be 
used as additional checks on the applicability of our model. 
The slope is matched quite well, even though this is not en- 
tirely trivial because it requires that the vertical structure 
that develops under the action of P B s for the self-gravitating, 
star-forming sheet of gas leads to star formation rates per 
unit area which are compatible with the observed slope of 
the Schmidt-law. Interestingly, the cut-off induced by the 



Figure 3. Star formation rate per unit area versus gas surface 
density in a self-consistent simulation of a disk galaxy which quies- 
cently forms stars. The symbols show azimuthally averaged mea- 
surements obtained for our fiducial choice of = 2.1 Gyr. The 
dashed inclined line gives the Kennicutt law of equation (25), and 
the vertical line marks the observed cut-off of star formation. 



best-fit value of to also lies approximately in the right loca- 
tion. It is presently unclear whether this has any profound 
significance, or whether it is just a fortunate coincidence 
in the present simple model. Recall that the cut-off in the 
model is induced by an imposed physical density threshold 
Pth for the onset of cloud formation, and that this density is 
tied to the value for the star formation timescale. 

Finally, we examine how well full three-dimensional sim- 
ulations of spiral galaxies obey the Kennicutt law that we 
used to set the star formation timescale. In Figure 3, we 
show azimuthally averaged measurements obtained for our 
fiducial choice of to =2.1 Gyr in a compound galaxy consist- 
ing of a dark halo, and a star-forming gaseous disk. There 
is good agreement with the corresponding analytic curve 
in Fig. 2, validating the numerical implementation of the 
multi-phase model in our simulation code. 



4 WINDS AND STARBURSTS 
4.1 Winds 

As summarised above, our multi-phase model leads to the 
establishment of a physically motivated and numerically well 
controlled regulation cycle for star formation in gas that has 
cooled and collapsed to high baryonic overdensities. Gas con- 
tained in dark matter haloes can thus cool and settle into 
rotationally supported disks where the baryons are gradually 
converted into stars, at a rate consistent with observations 
of local disk galaxies. In this model, the thickness and the 
star formation rates of gaseous disks are regulated by su- 
pernova feedback, which essentially provides finite pressure 
support to the star forming ISM, thereby preventing it from 
collapsing gravitationally to exceedingly high densities, and 
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also allowing gaseous disks to form that are reasonably sta- 
ble against axisymmetric perturbations. 

However, it is clear that the model we have outlined 
so far will not be able to account for the rich phenomenol- 
ogy associated with starbursts and galactic outflows, which 
are observed at both low (e.g. Heckman et al., 1995; Bland- 
Hawthorn, 1995; Lehnert & Heckman, 1996; Dahlem et al., 
1997; Heckman et al., 2000) and high redshifts (e.g. Frye 
et al., 2002; Pettini et al., 2000, 2001). This is because 
our multi-phase model by itself offers no obvious route for 
baryons to climb out of galactic potential wells after hav- 
ing collapsed into them due to cooling. Note that for the 
hybrid model of quiescent star formation we explicitly as- 
sume that the cold clouds and the hot surrounding medium 
remain tightly coupled at all times. The high entropy gas 
of supernova remnants is thus trapped in potential wells 
by being tied into a rapid cycle of cloud formation and 
evaporation. In principal, tidal stripping of enriched gas in 
galaxy interactions (Gnedin & Ostriker, 1997) could lead 
to a transport of enriched gas back into the low-density 
IGM. However, high-resolution simulations of galaxy colli- 
sions (Barnes, 1988; Barnes & Hernquist, 1992; Hernquist, 
1992, 1993) have shown that such dynamical removal of gas 
from the inner regions of galaxies appears to be rather in- 
efficient, especially for the deep potential wells expected for 
haloes in CDM universes (Springel & White, 1999) . 

On the other hand, it is becoming increasingly clear 
that galactic winds and outflows may play a crucial role 
not only in chemically enriching and possibly heating the 
IGM (Nath & Trentham, 1997; Madau et al., 2001; Aguirre 
et al., 2001a, b,c), in polluting the IGM with dust (Aguirre, 
1999a, b), and in enriching the intracluster and intragroup 
medium, but may also be an important mechanism in regu- 
lating star formation on galactic scales (Scannapieco et al., 
2000; Scannapieco & Broadhurst, 2001a,b). Since winds can 
reheat and transport collapsed material from the center of 
a galaxy back to its extended dark matter halo and even 
beyond, they can help to reduce the overall cosmic star for- 
mation rate to a level consistent with observational con- 
straints. Because radiative cooling is very efficient at high 
redshifts and in small haloes (White & Rees, 1978; White 
& Frenk, 1991), numerical simulations of galaxy formation 
typically either overproduce stars compared with the lumi- 
nosity density of the universe, or harbour too much cold gas 
in galaxies. The self-regulated model we present above will 
also suffer from this problem, because it does not drasti- 
cally alter the total amount of gas that cools. It is plausible, 
however, that galactic winds may solve this "overcooling" 
problem, provided that they can expel sufficient quantities 
of gas from the centres of low-mass galaxies. Removal of such 
low-angular momentum material may also help to resolve 
the problem of disk sizes being too small in CDM theories 
(Navarro & White, 1994b; Navarro & Steinmetz, 2000; Bin- 
ney et al., 2001). Note that semi-analytic models of galaxy 
formation must also invoke feedback processes that reheat 
cold gas and return to the extended galactic halo or eject it 
altogether. 

We are thus motivated to extend our feedback model to 
account for galactic winds driven by star formation. Winds 
have been investigated in a number of theoretical stud- 
ies (Mac Low & Ferrara, 1999; Aguirre et al., 2001a,b,c; 
Efstathiou, 2000; Madau et al., 2001; Scannapieco et al., 



2001, among others), but the mechanism by which galactic 
outflows originate is not yet well understood. In the star- 
forming multi-phase medium, it is plausible that not all of 
the hot gas in supernova remnants will remain confined to 
the disk by being quickly used up to evaporate cold clouds. 
Instead, supernova bubbles close to the surface of a star- 
forming region may break out of a disk and vent their hot 
gas into galactic haloes. As a result, a galactic-scale wind 
associated with star formation may develop. Note that this 
process does not necessarily require a prominent starburst, 
but could be a common phenomenon even with quiescent 
star formation (Efstathiou, 2000). In the latter case, winds 
may often not be strong enough to escape from all but the 
smallest galaxies, but they can nevertheless represent a cru- 
cial process for transporting metals and energy released by 
supernovae into the extended galactic halo. 

Unfortunately, our understanding of the small-scale hy- 
drodynamics that is actually responsible for driving winds 
is crude, although high-resolution simulations of starbursts 
(Mac Low & Ferrara, 1999) provide some clues about the 
physics. Here, we will not attempt to formulate a detailed 
theory for the production of winds, but will rather use a phe- 
nomenological approach that can easily be combined with 
our model for star formation and feedback. To this end, 
we will parameterise the winds in analogy to Aguirre et al. 
(2001b), who studied wind propagation and metal enrich- 
ment of the IGM using detailed analytic models applied to 
a sequence of simulation time slices. Their phenomenolog- 
ical model of winds was constructed to be physically and 
energetically plausible, to be simple, and to be constrained 
by observational data to the extent possible, an approach 
we will also follow. 

We start by assuming that the disk mass-loss rate that 
goes into a wind, M w , is proportional to the star formation 
rate itself, viz. 

Mw = vM+ , (26) 

where 7/ is a coefficient of order unity. In fact, Martin (1999) 
presents observational evidence that the disk mass-loss rates 
of local galaxies are of order the star formation rates them- 
selves, with M w /M* ~ 1 — 5, and without evidence for any 
dependence on galactic rotation speed. Note that M w just 
describes the rate at which gas is lost from a disk and fed 
into a wind. Whether or not this material will be able to es- 
cape from the halo will then depend on a number of factors: 
the velocity to which the gas is accelerated, the amount of 
intervening and entrained gas, and the depth of the potential 
well of the halo. If the wind is slow and the halo is of suffi- 
cient mass, the gas ejected from the disk will remain bound 
to the halo, and may later fall back to the star-forming re- 
gion, giving rise to a galactic "fountain." On the other hand, 
even a slow wind may escape dwarf galaxies and, at suffi- 
ciently high redshift, can potentially pollute a large volume 
without overly perturbing the thermal structure of the low- 
density IGM. Note that for values of r\ above unity, a wind is 
expected to greatly suppress star formation in galaxies that 
are of low enough mass to allow the wind to escape. We will 
adopt r\ — 2 in our fiducial wind model, consistent with the 
observations of Martin (1999). 

We further assume that the wind carries a fixed fraction 
X of the supernova energy. Equating the kinetic energy in 
the wind with the energy input by supernovae, 
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^Mw«w = X^snM*, (27) 
we obtain the wind's velocity when it leaves the disk as 



(28) 



We will treat \ as a further parameter of the wind model. 
While it cannot be larger than 1 in the present formalism, 
it is nevertheless to be expected that the wind energy con- 
stitutes a sizable fraction of the overall supernova energy 
input. We will assume \ = 0.25 in most of the test simula- 
tions discussed in this paper. Note that it should in principle 
be possible to use observations of the cosmic star formation 
history or the metal enrichment of the IGM to constrain 
this parameter. We do not explicitly model a distribution of 
wind velocities in this study. Of course, some variation of 
the local wind speed will be automatically produced by the 
dynamics of local gas interactions. 



the effective pressure model discussed so far shows a pos- 
sibility for run-away star formation, for purely dynamical 
reasons. Recall that the ISM is stabilised against gravi- 
tational collapse by the effective pressure provided in the 
multi-phase model. For sufficiently high densities, the cor- 
responding equation of state eventually becomes soft. More 
specifically, there is a certain overdensity, where the local 
polytropic index falls below a slope of 4/3. It is well known 
that barotropic gas spheres with an index below this value 
are unstable. We thus expect that once we assemble a suffi- 
ciently large amount of cold gas, the effective pressure will 
be insufficient to stabilise the cloud against further gravita- 
tional collapse to much higher densities, leading to a run- 
away conversion of the gas into stars. 



4.2 Starbursts 

We have shown that our hybrid model leads to an efficient 
self-regulation of star formation, which can be understood in 
terms of an effective equation of state. In general, the model 
is thus best viewed as describing quiescent star formation, 
where the star formation rate will in general only gradually 
accelerate when the gas density is increased, in accordance 
with the empirical Schmidt law. The winds we have intro- 
duced as an extension of this model will not change this 
picture. We specifically postulate that a wind leaves a star- 
forming region without significantly perturbing it dynami- 
cally, and our numerical implementation of wind formation 
is designed to ensure this behaviour. Winds thus merely re- 
duce the amount of gas available for star formation. De- 
pending on the depth of the potential well, the gas may or 
may not come back to the galactic disk at a later time and 
become available for star formation again. 

It is interesting to note, however, that we expect the qui- 
escent mode of star formation to eventually become "explo- 
sive" (i.e. very rapid) for sufficiently high gas densities. Phys- 
ically, it is plausible that self-regulation should break down 
at high densities. In the self-regulated regime, cold clouds are 
constantly being formed and evaporated. If clouds are not 
replenished by cooling, they will be consumed on a timescale 
t c = U/(J3A) = (p/pth) 3/W to/(l3Ao). This timescale t c de- 
scribes the rate at which clouds are reprocessed. At the on- 
set of star formation, it is about a factor 100 shorter than 
the star formation timescale itself, on the order of a few 
times 10 7 years. However, at higher densities, clouds are re- 
processed more slowly, and eventually t c will become larger 
than the maximum lifetime of individual clouds, which is 
estimated to be as high as 10 8 yr for giant molecular clouds. 
At this point, it will become difficult to maintain tight self- 
regulation because the clouds will not survive long enough 
to await evaporation. Instead they may all collapse and form 
stars on the timescale of their lifetime, i.e. the timescale of 
star formation will suddenly become shorter in this regime 
and deviate from the scaling we have assumed so far. 

Since it is unclear how this transition to accelerated 
star formation proceeds in detail, we refrain from modelling 
it explicitly in this study. However, we remark that already 



5 NUMERICAL IMPLEMENTATION 

The hybrid multi-phase model outlined in Section 2 can be 
represented numerically in a number of different ways. Yepes 
et al. (1997) used an Eulerian hydro code, where the prop- 
erties of the flow are discretised on a mesh. The mesh size 
provides a natural scale for the coarse- graining procedure 
inherent in our model. One may also discretise the mass, 
as is done in the Lagrangian SPH technique. Huffman & 
Pharasyn (1999) presented a first implementation of the ap- 
proach of Yepes et al. (1997) using SPH. 

We also employ SPH to implement our model, in the 
form of a modified and extended version of the parallel 
TreeSPH code GADGET (Springel et al., 2001b). Through- 
out, we will use the novel formulation of SPH that was re- 
cently derived by Springel & Hernquist (2002a), which man- 
ifestly conserves energy and entropy when appropriate, and 
greatly reduces numerical inaccuracies that otherwise can 
arise in poorly resolved cooling flows. 

On a technical level, we have actually implemented two 
different methods for representing our star formation model 
in the code. In an "explicit" treatment, which we will de- 
scribe first, each SPH particle can have three different mass 
components, describing "ordinary" ambient gas, cold clouds, 
and collisionless stars. The various processes of cloud for- 
mation, cloud evaporation and star formation are then ex- 
plicitly followed in terms of mass exchange between these 
components of the hybrid particles. 

However, because the self-regulation timescale is typi- 
cally short compared to the dynamical and star formation 
timescales, the mass fraction contained in cold clouds at a 
given density can be obtained with good accuracy by just 
assuming the equilibrium value expected for self-regulation, 
where the latter can be worked out analytically. We can 
thus replace the explicit treatment with a simplified method 
based on these equilibrium values. Star particles can then be 
formed stochastically from the reservoir of clouds, thereby 
avoiding an artificial coupling of gas and collisionless stellar 
material. This idea forms the basis of our second method for 
implementing the multi-phase model, to be described after 
the explicit treatment. 
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5.1 Explicit implementation of the multi-phase 
model 

In our explicit variant for implementing the multi-phase 
model, we describe each SPH hybrid-particle by a total mass 
m, a mass m* in stars, a mass m c in cold clouds, and an 
internal energy Uh in the ambient gas. The mass in the "or- 
dinary" gas phase is therefore given by mj, — m — m* — m c . 
Note that we here build up stellar mass gradually in each 
star-forming particle, unlike in the effective treatment dis- 
cussed in the next section, which avoids the implied tempo- 
rary coupling of gas and collisionless stellar material. Also 
note that most SPH particles in a cosmological simulation 
will never reach the high densities required for star forma- 
tion, and will therefore have m = mu, with Uh describing 
their normal thermal reservoir. 

We compute the hydrodynamical forces in the usual 
manner, but take the pressure to be 



= m* + (1 — P)Am s f 



(37) 



P = (7 - l)p(m h u u + m c u c )/(m h + m c ), 



(29) 



and compute kernel-estimated densities with the gas mass 
m — m* of each particle only. The stellar masses m* add 
inertia to the particles and contribute to the source func- 
tion of the gravitational field. Star formation is then treated 
in the code on a particle by particle basis in the following 
manner: 

(i) For each timestep At, we first work out the mass of 
clouds that turns into stars and produces supernovae: 

At 



Am,f = m c 



(30) 



(ii) If the thermal instability is active (/ = 0), we then 
compute the new mass in the hot phase. To this end, we 
consider the Lagrangian analogue of equation (6): 



Am h = /3(A+ l)Am Bf - -^-^At, 



(31) 



Uh - U c Ph 

and we solve this equation implicitly for the new mass m' h 
in the following way: 

' , lU A net (p h m' h /m h ,u h ) m h 
m h — mh + p{A + 1) Am s f — At(32) 

Uh - u c p h 

This guarantees stability even when rrih changes rapidly. In 
order to improve the robustness of the integration, we how- 
ever do not allow m' h to fall below half the value of mj, in a 
single step. If the thermal instability is not active, we simply 
have m' h = m h + fi(A + l)Am s f . 

(iii) Once m' h is determined, we can compute the mass 
that is actually transferred from the hot to the cold phase 
due to the thermal instability alone: 

Amu = m h + (3 (A + l)Am sf - m' h . (33) 

Note that Amxi = if the thermal instability is not oper- 
ating. The mass evaporated by supernovae is 



Am. 



(3AAm a( . 



(34) 



(iv) The new masses of cold and hot phase can now be 
obtained as: 



m' c = m c — Am sf + Am T i — Am evap 
m'h = m h + PAm s { - Am T i + Am ovap 
And the new stellar mass is: 



(35) 
(36) 



Note that the total mass is conserved, m' h + m' c + m'+ — 
mh + rn c + rn+, by construction. 

(v) If the thermal instability is operating, the rate of 
change of internal energy of the ambient medium is given 
by the normal adiabatic term plus the contribution 

Au h _ /3Am s{ [m sn + (A + l)(u c - u h )] 
At 



(38) 



m h At 

from star formation and feedback. Otherwise, if ordinary 
cooling is active, the rate of change of internal energy ad- 
ditionally gets a contribution from the normal cooling func- 
tion, with the latter being solved implicitly for the new tem- 
perature at the end of the timestep (Springel et al., 2001b). 

Note that the coupling of the three phases can cause 
difficulties in certain situations. For example, if a hybrid 
particle is dynamically stripped from a star forming region 
and moves to an environment of lower density, the particle 
will at first have most of its gas mass still bound in cold 
clouds. These clouds will quickly be evaporated by residual 
star formation of the particle, provided its density has not 
fallen so rapidly that the star formation timescale has be- 
come very long. If this should happen, cloud material may 
escape to very low density, and survive there for a long time. 
In order to prevent this, we assume that all clouds are evap- 
orated instantly if the density of a particle falls below p th , 
with the required energy extracted from the reservoir Uh- 
Star formation will thus be strictly confined to the regions 
p > pth, and all gas of lower density is just "ordinary" gas. 

Somewhat more problematic is the tight coupling of the 
stellar material that has formed to the mass that still resides 
in gas. This is in principle unphysical, even though it may 
be a reasonable approximation in many cases, because the 
stars and the star-forming medium tend to move together, 
except in highly dynamic situations, as, for example, in ma- 
jor mergers. Nevertheless, it is clear that the explicit treat- 
ment described above requires a mechanism that eventually 
decouples the formed stellar material from the gas, other- 
wise one could obtain at late times large numbers of "gas" 
particles whose masses are dominated by their stellar com- 
ponents. In order to achieve such a decoupling, we turn a 
hybrid particle whose stellar mass has grown above its gas 
mass, i.e. for m* > 0.5m, into a star particle of mass m*. 
The rest of the mass, consisting of clouds and ambient gas, 
is distributed kernel-weighted among the SPH neighbours 
of the particle that is being converted. We distribute clouds 
and ambient mass separately, as well as the thermal energy 
in the ambient medium. 

This mechanism maintains a roughly constant mass res- 
olution, and it keeps the total particle number exactly con- 
stant. However, it is clear that subtle dynamical effects due 
to the temporary coupling of gaseous and stellar material 
may remain. This is one of the reasons why we prefer a sim- 
plified, "effective" implementation of the multi-phase model 
which we describe next. 



5.2 Simplified treatment assuming self-regulation 

According to equation (10), the ambient hot medium evolves 
towards an equilibrium temperature on a timescale given 
by equation (12), provided that the density is sufficiently 
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high for star formation to occur. This time is short com- 
pared to the star formation timescale; i.e. conditions for 
self-regulation are achieved quickly. In order to simplify the 
computation of the gasdynamics, it should, therefore, be a 
good approximation to assume that the conditions of self- 
regulated star formation always hold. In this case, the tem- 
perature of the ambient hot medium is simply given by equa- 
tion (11), and the mass fraction in clouds follows from equa- 
tions (16) and (18). This then also determines the star for- 
mation rate and the effective pressure (19) of each particle. 
In this way, we avoid treating the mass exchange processes 
explicitly, simplifying the code, making it faster, and reduc- 
ing its storage requirements. 

Additional advantages follow if the build-up of the stel- 
lar component is not described "smoothly" as in the ex- 
plicit approach, but instead probabilistically with expecta- 
tion value consistent with the star formation rate. The star 
formation rate of an SPH particle of current mass in is given 
by Mi, = (1 — P) x m/ti,. Given a timestep At, we spawn a 
new star particle of mass m* = mo /N g once a random num- 
ber drawn uniformly from the interval [0, 1] falls below 



m 
m* 



exp 



(l-/3)xAt 



(39) 



The initial phase-space variables of the new particle are 
copied from the gas particle, which in turn is reduced in 
mass by m*. Here, mo is the initial gas mass of each SPH 
particle, and N g is an integer which gives the number of 
"generations" of stars each gas particle may form. If a gas 
particle has already spawned N g — 1 stars, it will simply 
be turned into a star particle should it become the site of 
another star formation event. 

Note that all the star particles will have identical masses 
in this approach. For choices of N g in the range between 
1 and about 4, one achieves good mass resolution for the 
stellar component, while simultaneously maintaining a rea- 
sonably constant mass resolution in the gaseous phase. For 
N g > 1 the total number of particles increases due to star 
formation, but there is a firm upper bound to this number, 
and if only a small fraction of the baryons are turned into 
stars, as in cosmological simulations, the relative increase is 
moderate. 

More important, there is no artificial dynamical cou- 
pling between the gas particles and the stars in this ap- 
proach, because all the stellar mass is contained in indepen- 
dent star particles at all times, unlike in the explicit treat- 
ment, where each component is tied to hybrid gas-star parti- 
cles. Another advantage is that the distribution of formation 
times of the star particles directly reflects the underlying 
star formation history, without any bias, thereby simplifying 
the analysis and interpretation of simulation outputs. The 
equal masses of star particles also keep possible mass seg- 
regation effects due to collisional relaxation to a minimum, 
and in the limit of very large N g one naturally approaches 
the limit of continuous star formation. 



AM* locked up in long-lived stars, the mass in metals re- 
turned is AM mct = pAM*, where p is the yield. Provided 
that these metals are well mixed with the local gas of mass 
M g , and provided no gas enters or leaves the system (i.e. a 
closed-box) the metallicity Z = Mm et /M g will increase by 
AZ = pAM+ / Mg due to the formation of the stars. This 
already takes into account that some of the metals in the 
gas are also lost to the forming stars. 

Here, we assume that each gas element locally behaves 
as a closed box, with metals being instantaneously mixed 
between clouds and ambient gas. The metallicity of a patch 
of the ISM will therefore increase during one timestep by 



AZ = (l-/3)px 



At 



(40) 



where x is the local fraction of gas contained in cold clouds. 
For gas below the threshold density for star formation, we 
have x = and the metallicity remains constant. 

It is straightforward to follow the growth of the metal- 
licity of each gas particle using equation (40), both for the 
explicit treatment of the multi-phase structure as described 
in section 5.1, or for the simplified "effective" method de- 
scribed in section 5.2. However, depending on the approach 
used, it is more problematic to define the metallicity of the 
star particles. In the explicit treatment, the stellar compo- 
nent is built up gradually, with stars being formed from gas 
of varying metallicity. The independent star particles that 
are eventually produced will then have a mass-weighted av- 
erage of the metallicities of the gas elements from which 
they were formed. This complicates the interpretation of the 
metallicity distribution of these star particles, because it will 
not directly reflect the metallicity distribution obtained if 
star formation were instead always to immediately generate 
independent star particles. 

In the probabilistic method of generating star parti- 
cles, this problem does not arise. The metallicity of a newly 
spawned star particle is just given by the current metallicity 
of the star-forming gas particle. In this method, the expec- 
tation value of the total mass in metals found in gas and 
stars at any given time will be equal to the yield times the 
total mass in stars. Note that the effective model also avoids 
any metal diffusion, which can occur at a small level in the 
explicit treatment where remaining gas and metals are dis- 
tributed among neighbouring particles once a gas particle 
is dissolved, with its hybrid stellar component being turned 
into a collisionless particle. 

Although there is no efficient mechanism for metal diffu- 
sion in our model, metals can of course be transported along 
with the gas flow. In particular, winds may deposit metals 
far from where they were produced, leading, for example, to 
enrichment of the haloes of galaxies, or of the IGM. How- 
ever, currently we do not yet include metal-line cooling in 
our code, hence the metallicity acts fundamentally only as a 
tracer variable without having any dynamical consequences. 
In future work, we plan to investigate effects that arise when 
this part of the model is refined. 



5.3 Metal enrichment 

Stellar winds and supernova explosions of massive stars en- 
rich surrounding gas with metals (i.e. all elements heav- 
ier than helium). We assume that for each mass element 



5.4 Wind formation 

As discussed in Section 4, we assume that star formation can 
be accompanied by mass-loss from star-forming regions, giv- 
ing rise to winds. There is strong observational evidence that 
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such winds exist, but the detailed physics underlying this 
phenomenon is unclear. But, even if there were a firm theo- 
retical basis for understanding these winds, we currently lack 
the spatial resolution in our simulations to directly model 
the relevant physical processes. We are thus led to adopt a 
phenomenological treatment of how such winds are gener- 
ated. At a technical level, we proceed in a similar way to 
the probabilistic method for spawning star particles. 

During a timestep At, a gas particle is added to the 
wind if a uniformly distributed random number out of the 
interval [0, 1] falls below 



exp 



77(1 -13) xAt 



(41) 



If this is the case, we modify the particle's velocity v ac- 
cording to 



(42) 



where v w is the wind velocity given by equation (28). For 
the unit vector n, we either select a random direction on the 
unit sphere (isotropic wind), or we choose a random orienta- 
tion along the direction v x V0, where <f> is the gravitational 
potential ( "axial" wind) . If the latter scheme is chosen, the 
wind particles are preferentially ejected along the rotation 
axis of a spinning object. This mimics the expected prefer- 
ence of wind ejection to occur in a direction orthogonal to 
rotationally flattened, disk-like galaxies. However, even for 
an isotropic wind the presence of a dense disk will lead to 
a bipolar wind pattern, simply because the wind can prop- 
agate much more easily in the direction orthogonal to the 
disk. Note that because of the random orientation of the mo- 
mentum kick along the direction n, the expectation values 
of momentum and angular momentum remain unchanged. 
While the momentum is not strictly balanced when a single 
particle is put into the wind, it is hence statistically con- 
served. Also, the expected average increase of the specific 
kinetic energy of a particle that enters the wind is given by 
Vto/2, as desired. 

An interesting problem arises because of the finite thick- 
ness of a star-forming region. A real wind plausibly origi- 
nates from a region close to the surface of a star-forming 
disk, allowing it to leave without greatly impacting ongoing 
star formation. In the simple model developed here, we do 
not try to restrict wind formation to a surface layer. Instead, 
we allow all the SPH particles in the star-forming region to 
enter the wind if chosen by the probabilistic criterion. How- 
ever, wind particles from inner parts of star-forming disks 
would normally be stopped again by other particles in the 
dense disk, resulting in rapid dissipation and thermalisa- 
tion of the kinetic energy. In addition, the momentum input 
would lead to a strong perturbation and, in extreme cases, 
to a disruption of the disk. This is in disagreement with 
our assumption that the wind can escape from the dense, 
star-forming phase without directly affecting it. Only out- 
side the disk, would gasdynamical interactions within the 
halo stop the wind. To allow such behaviour, we "decouple" 
a spawned wind particle for a brief time from hydrodynamic 
interactions, i.e. it neither exerts nor receives hydrodynamic 
forces during this period. The particle is however included 
in the computation of the gas density, and for gravity. The 
full hydrodynamic interactions are enabled again once the 
density of the particle has fallen to 0.1 pth, or once a time of 
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Figure 4. Star formation rates in simulations of disks growing 
within haloes of different masses, with or without the inclusion of 
a wind. From top to bottom, the panels show results for haloes 
of differing total mass: M vir = 10 10 /i _1 M Q , 10 11 7i _1 M , and 
10 12 /i _1 M0, respectively. The inclusion of winds always leads to 
a net suppression of star formation and hence to a lower curve, 
but the importance of this effect is a strong function of halo mass. 
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At = 50 Myr has elapsed, whichever happens earlier. This 
allows the wind particle to travel "freely" up to a few kpc 
until it has left the dense star-forming phase. In this way, 
we can mimic quite well the strong, but quiescent mass- 
loss from a star- forming region. We have checked that the 
strength of the wind and the efficiency of wind escape are 
insensitive to the detailed values chosen for the parameters 
of this decoupling prescription. 



6 TESTS AND RESULTS 

The feedback model discussed above has widespread impli- 
cations for cosmological simulations, influencing nearly all 
predicted properties of galaxy populations. However, it is 
beyond the scope of this paper to fully investigate all these 
consequences in detail. Rather, we focus here on the method- 
ology of the modelling technique. The simulations we present 
below were selected as test cases to allow us to validate our 
approach, and to give us insight into key effects of the feed- 
back model in a cosmological context. 

More specifically, we focus on two small sets of simula- 
tions, one designed to study the problem of disk formation in 
a highly idealised situation, and the other selected as a full- 
blown, yet rather small cosmological simulation of structure 
formation. Results of modelling with higher dynamic range 
will be presented in due course. 



6.1 Isolated spiral galaxies 

Currently, simulations of cosmic structure formation still fail 
to match observed properties of spiral galaxies, although 
some progress has been achieved recently using improved 
feedback prescriptions (e.g. Sommer-Larsen et al., 2002). 
One cause of the difficulty has been that collapsed gas loses 
too much angular momentum to the dark matter, so that 
disks are too compact when compared with observations. 
On the other hand, the spin distribution of dark matter 
haloes is well understood both analytically and numeri- 
cally. Simple analytic models indicate that disk sizes in 
such haloes should match the observed distribution, pro- 
vided that baryons conserve most of their specific angular 
momentum during collapse (Mo et al., 1998; van den Bosch 
et al., 2001; van den Bosch, 2001; van den Bosch et al., 2002). 

In our tests here, we will not be concerned with these 
issues directly. Instead, we examine the inside-out formation 
of a disk galaxy under idealised settings which allows us to 
cleanly study the properties of our self-regulated model for 
star formation, with or without winds. To this end, we con- 
sider a series of simulations which effectively follow the Fall 
& Efstathiou (1980) picture of the formation of a disk galaxy. 
We set up an NFW-halo (Navarro et al., 1996, 1997) in iso- 
lation, with the dark matter and gas initially in virial equi- 
librium. When evolved in time using adiabatic gas physics, 
our initial conditions are nearly stable, showing only negli- 
gible secular evolution over a Hubble time. However, when 
the gas is allowed to cool radiatively, the gas in the center of 
the halo quickly loses its pressure support and settles into a 
rotationally supported disk that grows from inside out, with 
a size that depends on the initial angular momentum of the 
gas. Note that due to the axisymmetric set-up we have cho- 



sen, angular momentum transport between the gas and the 
dark matter is essentially absent in these simulations. 

For definiteness, we have chosen a halo of mass 
M v i r = 10 12 ft _1 M 01 with 10% of the mass being in bary- 
onic form. The initial angular momentum J of the halo 
can be described in terms of the spin parameter A = 
J\E\ 1/2 /(GM% 2 ), for which we adopt a value A = 0.1 to 
produce a large disk. The relative distribution of angular mo- 
mentum within the halo was chosen under the assumption 
that the specific angular momenta j(r) of spherical shells are 
all aligned, and that their magnitude is given by the fitting 
formula 

of Bullock et al. (2001). Within each shell, we distribute the 
angular momentum as if the shell were in solid body rota- 
tion. For s = 1, the choice we adopt here, this implies that in 
the xy-pl&ne, the initial azimuthal streaming velocity is pro- 
portional to the circular velocity squared. We have also sim- 
ulated haloes with masses of 10 11 /i _1 M , and 10 10 /i _1 M©. 
They are simply scaled down versions of the 10 12 /i _1 M 
halo, differing by factors of 10 in mass, and by 10 1/3 ~ 2.15 
in length and velocity scale. 

In order to simplify these test simulations further, we 
have represented the dark matter halo in most of our runs by 
a static gravitational potential. This neglects the contraction 
of the halo due to the infall of the gas, but for the purposes of 
the present analysis this effect is not important. Simulations 
that include a live dark halo give very similar results, but 
are more expensive, especially when one wants to reduce the 
additional particle noise from the dark halo to a negligible 
level. 

For each of the three halo masses, we run simulations 
with and without a wind. Typically we used 40000 self- 
gravitating SPH particles in each case, except for a number 
of resolution tests which we will discuss separately below. 
The parameters of our multi-phase model were set to our 
set of fiducial values, i.e. to = 2.1 Gyr, and for the wind 
sector, f] — 2 and x = 0.25. 

In Figure 4, we show the star formation rates in these 
simulations as a function of time. In all cases, the rapid ini- 
tial development of the gas disk is accompanied by an abrupt 
rise of the star formation rate to a maximum soon after the 
start of the simulations. Thereafter, the star formation rate 
declines nearly exponentially for some time, but the decline 
becomes increasingly slower at later times. In general, the 
models with winds yield a lower star formation rate at all 
times, but the magnitude of the suppression is a strong func- 
tion of halo mass. The 10 12 ft"'M0 halo forms stars nearly 
unaffected by a wind, but the 10 h~ Mq halo and particu- 
larly the 10 10 /i _1 Mq halo suffer a strong reduction in their 
star formation rates. 

The cause for these differences between the effects of 
the wind in haloes of different mass becomes readily appar- 
ent when one examines the flow of the gas in more detail. In 
Figures 5, 6, and 7, we show the time evolution of the veloc- 
ity field in the xz-pl&ne, which is a slice orthogonal to the 
gas disk, with the rotation axis along the z-axis. We indicate 
the gas velocity field with arrows, overlaid on a colour-scale 
map showing the gas density in the slice. 

For the massive 10 12 h~ 1 M & halo, the wind is never able 
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Figure 5. Time evolution of the gas flow in a halo of total mass M v ; r = 10 Ii _1 Mq. The velocity field is represented by arrows and 
the logarithm of the gas density is indicated as a colour-scale. Labels give the elapsed time in /i -1 Gyr since the start of the simulation. 
A wind of speed 242kms~ 1 is included in this model, but it is stopped already very close to the disk. 
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Figure 6. Time evolution of the gas flow in a halo of total mass M v ; r = 10 11 Ii -1 Mq. The velocity field is represented by arrows and 
the logarithm of the gas density is indicated as a colour-scale. Labels in each panel give the elapsed time in h~ 1 Gyr since the start of 
the simulation. A wind of speed 242kms — 1 is included in this model, somewhat slower than the escape speed of v eS c ~ 280 km s -1 from 
V^/U&lo^As, a^ijc^u^^^a^c^cagcsQijOj substantial heights in the halo, but falls back later onto the disk in a galactic "fountain". 
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Figure 7. Time evolution of the gas flow in a halo of total mass Af v i r = 10 10 /i _1 Mq. The velocity field is represented by arrows and 
the logarithm of the gas density is indicated as a colour-scale. Labels in each panel give the elapsed time in h _1 Gyr since the start of the 
simulation. A wind of speed 242 kms -1 is included in this model, considerably higher than the escape speed of v CBC ~ 130 kms -1 from 
this halo. As a result, a galactic super- wind develops which blows out of the galaxy, entraining a significant fraction of the gas from the 
halo. 
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Figure 8. The same as Figure 7, except that the simulation included isotropic winds instead of explicitly modelling an outflow that is 
preferentially oriented along the rotation axis of the disk. A bipolar outflow pattern orthogonal to the disk still develops for the isotropic 
wind due to the dense disk that forms in the xy-plane, and due to the non-isotropic inflow pattern. Overall, the gas flow is considerably 
less ordered in this case, but it is qualitatively still quite similar to the model with axial winds. 
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Figure 9. Stellar disk that formed after 3 Gyr in a halo of mass M v ] r = 10 12 /i _1 Mq, seen both face on and edge on. On the left, the 
simulation does not include winds, on the right it does. The colour-contours in the pictures have been set such that they contain roughly 
the same fraction of the total light in each case. 



to propagate very far from the star-forming disk. Instead, 
any ejected material is quickly stopped by the ram pressure 
of ambient gas, and, even more important, by the strength 
of the gravitational force field. Note that the wind speed of 
242 kms -1 is comparable to the circular velocity of this halo, 
but it is still much smaller than the halo's escape velocity. 

However, the evolution is already somewhat different 
in the 10 h~ M© halo, which is shown in Figure 6. In 
this case, the wind is still not energetic enough to unbind 
material from the halo, which has an escape velocity of 
v csc — 280kms _1 . However, ejected gas is able to rise to 
substantial heights in the halo and to displace infalling ma- 
terial, before it eventually falls back to the disk, forming 
a "galactic fountain" (Shapiro & Field, 1976). This process 
slows down the processing of gas into stars, and it mixes met- 
als into the halo, which thus becomes gradually enriched. 

Finally, if the potential well is shallow enough, the wind 
can escape from the galaxy entirely. This is the case for the 
10 10 /i _1 Mq halo, as shown in Figure 7. The wind veloc- 
ity is now substantially larger than the escape velocity of 
?; csc ~ 130kms _1 , allowing the wind to "blow out" into the 
intergalactic medium, where it deposits both metals and its 
residual kinetic energy. The wind can also entrain some of 
the halo's gas, thereby further reducing the infall rate onto 
the star-forming disk. Note that our wind model has been 
specifically constructed to be "quiescent", i.e. to leave the 
star-forming ISM intact. The wind does not blow away the 
ISM entirely, but it instead leads to quiescent mass-loss with 
strength determined by the star formation rate. This re- 



moval of baryons, together with the interception of infalling 
material, leads to a strong reduction of the star formation 
rate in low-mass haloes. 

In the three wind runs discussed above, "axial" wind 
formation has been used. In the small-mass galaxy shown 
in Figure 7, this leads to a highly collimated, bipolar out- 
flow. In Figure 8, we show the same model, but computed 
with isotropic winds. In this case, a bipolar outflow pattern 
orthogonal to the disk still develops, but due to the dense 
disk that forms in the xy-plane and the non-isotropic in- 
fall pattern of the gas. Overall, this results in a gas flow 
that is qualitatively still quite similar to the model with ax- 
ial winds, and the resulting star formation rate is also very 
similar. In practice, even in this extreme case of a thin cen- 
trifugally supported disk, it makes little difference whether 
isotropic winds or the axial model are used. 

It is also interesting to study the stellar disks that de- 
velop in these simulations. In Figure 9, we show the lumi- 
nosity density of the stellar disk that formed in the two 
10 12 /i _1 Mq simulations after a time of 3 Gyr. As expected, 
a thin stellar disk has grown. Interestingly, the disk morphol- 
ogy appears nearly unaffected by whether or not a wind is 
included. In both cases, the disks are quite smooth, as are 
the gaseous layers that form the stars. This is a result of the 
pressurisation of the gas in our multi-phase model. Note that 
if the gas was treated as a single-phase medium, it would es- 
sentially behave as an isothermal gas at temperature 10 4 K. 
The disk would then be rather unstable to gravitational per- 
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Figure 10. Radial surface brightness profiles for disks that formed in isolated haloes after 3Gyr. On the left, we show results for a 
halo of total mass 10 10 /i _1 M0, both for a run without winds (stars) and one with winds (diamonds). The dotted lines are exponential 
profiles, fitted under the constraint that they give the same total light as the entire simulated disk. The exponential scale lengths of the 
disks are 0.51 and 0.58 h~ x kpc, respectively. On the right, the same results are shown for a halo of total mass 10 12 /i _1 Mq, with the 
exponential fits giving scale lengths 2.1 and 1.9 h~ kpc, respectively. 




-2-1012 -2-1012 

x [ h~ l kpc ] x [ /j" 1 kpc ] 



Figure 11. The stellar disk that formed after 3 Gyr in a halo of mass M v j r = 10 10 h~ 1 Mq, seen both face on and edge on. On the left, 
the simulation does not include winds, on the right it does. The colour-contours in the pictures have been set such that they contain 
roughly the same fraction of the total light in each case. Note that the the morphology of the two disks appears similar, but the model 
that includes the wind has only 36% of the luminosity of the run that does not. 
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Figure 12. Convergence study of the star formation rate in four 
runs that differ in total by a factor 64 in mass resolution. The 
simulations include winds and describe the formation of a disk in 
a 10 10 /i _1 Mq halo. We show results for simulations with 40000 
(solid), 160000 (dashed), 10000 (dot-dashed), and 2500 (dotted) 
particles. The gravitational softening used in the simulations was 
varied in proportion to the cube root of the mass resolution, rang- 
ing from 0.03 to 0.12/i- 1 kpc. 

turbations and would be bound to quickly break up into 
many lumps. 

The radial surface brightness profiles of the disks are 
shown in Figure 10. In the main body of the galaxy, the 
profile is reasonably well fitted by an exponential, but there 
is a clear deviation towards higher luminosities in the centre. 
This is probably a consequence of the initial angular momen- 
tum distribution that we adopted. For similar assumptions 
about the distribution of initial angular momentum, ana- 
lytic computations of disk formation predict a similar effect 
(van den Bosch, 2001). 

Figure 10 also shows how the influence of the wind 
depends on halo mass. While the stellar profile for the 
10 12 /i _1 Mq halo is nearly unaffected by the winds, the small 
10 10 /» _1 M0 halo suffers a reduction in its surface brightness 
by nearly a factor of 3. However, the radial scale lengths are 
essentially unaffected, as are the disk morphologies. The lat- 
ter can be seen in more detail in Figure 11, where we show 
the projected stellar mass densities for the low-mass halo. 
In comparison to Fig. 9, it is clear that the ratio of disk 
scale height to radial scale length is larger for the smaller 
halo. This is a direct consequence of the effective equation 
of state of the multi-phase medium. The effective pressure 
stabilises the gas disks vertically against gravity, but in do- 
ing so it imposes a certain scale, i.e. the scale height depends 
only weakly on surface mass density and does not fall below 
~ 150 pc in our model. 

Finally, we briefly examine the sensitivity of our numer- 
ical results to the mass resolution employed. This is a par- 
ticularly important test, not only to assess numerical con- 
vergence in general, but also for establishing the suitability 
of the method for cosmological simulations of structure for- 
mation. In such calculations, haloes of different mass are 
invariably resolved with a different number of particles. If 



the numerical results for star formation and feedback de- 
pend sensitively on mass resolution, differential resolution 
effects would be introduced that could severely compromise 
an interpretation of the results. Note that this also implies 
that the free parameters of any feedback model that behaves 
well in this respect should be independent of resolution. 

By construction, we expect the model described in this 
paper to meet these requirements, primarily because the 
feedback and star formation schemes are formulated in terms 
of an effective subresolution model which does not require an 
explicit reference to the particle formalism. We test this con- 
tention by repeating the wind-simulation of the 10 10 h~ 1 M.Q 
halo at a number of different mass resolutions, both with 
smaller and larger particle numbers than the 40000 used 
above. In particular, we selected mass resolutions of 2500, 
10000, and 160000 particles, thereby obtaining a series of 
runs that stretches a dynamic range of 64 in mass reso- 
lution. We also varied the 3D spatial resolution (i.e. the 
gravitational softening length) in proportion to the cube 
root of the mass resolution, as it is typically done in cos- 
mological simulations of structure formation. In Figure 12, 
we compare the star formation rates measured as a function 
of time in these simulations. Reassuringly, the agreement 
is rather good between all four runs, confirming that the 
model outlined above is numerically well posed. The for- 
mulation in terms of an effective sub-resolution model has 
a well-specified continuum limit, and we do not expect the 
results of this test to change even when the resolution is in- 
creased indefinitely. This means, however, that much higher 
resolution would not begin to eventually resolve the "true" 
spatial structure of the ISM - in order to do this faithfully, 
the effective model would have to be dropped. 

6.2 Cosmic star formation history 

We now employ our multi-phase model in full cosmolog- 
ical simulations of structure formation. In this paper, we 
restrict ourselves to highlighting some of the main conse- 
quences of the model with respect to the overall star forma- 
tion rate and the metal enrichment of the IGM. For these 
purposes, it is sufficient to consider moderate-sized simu- 
lations of a relatively small volume. In particular, we work 
with 2 x 50 3 dark matter and SPH particles in a box of length 
11.3 /i _1 Mpc per side. This gives a mass resolution of m gas = 
1.28 x 10 8 /i _1 M in the gas and m dm = 8.33 x 10 8 ft _1 M 
in the dark matter, for a canonical ACDM cosmology with 
parameters Q = 0.3, s2a = 0.7, D. b = 0.04, h = 0.67, and 
as = 0.9. We have run several of these simulations for dif- 
ferent parameter choices, and in addition one run where we 
stepped up the mass resolution by a factor of eight to 2 x 100 3 
particles. 

Note that the simulation volume in these runs is much 
too small to be representative of the universe as a whole, 
and in fact, the largest modes in the box become non-linear 
by z — 0. Nevertheless, the simulations suffice to identify the 
systematic effects that arise from our feedback schemes, and 
from the winds, in particular. However, it should be kept in 
mind that simulations of larger size with comparable or bet- 
ter mass resolution are required to arrive at quantitatively 
accurate results for cosmological expectation values such as 
the mean luminosity density. 

In Figure 13, we show the cosmic star formation history 
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Figure 13. Comoving star formation history of two small cosmo- 
logical simulations, using the hybrid multi-phase model without 
winds for star formation timescales = 1 Gyr and t ( * = 4.2 Gyr, 
respectively. For a smaller value of tj, the peak of the star for- 
mation rate moves to higher rcdshift, and the decline towards 
the present epoch proceeds faster. However, the integrated mass 
of stars that forms up to z = is insensitive to ij: the short- 
timescale model turns 18.2% of the baryons into stars, the long- 
timescalc model 18.4%. 



Figure 14. Evolution of the cosmic star formation density for 
different simulations. The solid line gives the star formation rate 
per comoving volume for our fiducial 2 X 50 3 simulation, with 
parameters = 2.1 Gyr, rj = 2, and \ = 0.25. The dot-dashed 
line is for the same model, but at a higher resolution of 2 X 100 3 . 
The dashed line shows a model with a much more energetic wind, 
specified by rj = 0.5 and x = 1> while the dotted line shows the 
result for a run without winds. 



for two of our wind-less 2 x 50 3 simulations, differing only in 
the value adopted for the star formation timescale tg. While 
a larger value of this parameter shifts star formation towards 
later times, the integrated star formation history is insen- 
sitive to the value of tj. Ultimately, if only the quiescent 
model of star formation is considered, the mass of stars that 
forms is approximately given by the amount of baryons that 
can cool, because the gas consumption timescale is signifi- 
cantly shorter than the age of the universe. It is, however, 
well known that cooling alone is quite efficient (e.g. White 
& Rees, 1978; White & Frenk, 1991; Katz & White, 1993; 
Pearce et al., 1999; Lewis et al., 2000; Dave et al., 2001; 
Balogh et al., 2001; Springel & Hernquist, 2002a), leading 
to a collapse fraction well in excess of observational bounds 
on the amount of baryons in stars and cold gas. For exam- 
ple, Fukugita et al. (1998) use optical data to constrain the 
fraction of baryons in the form of cold gas or stars to be in 
the range of 6.2% to 16.7%. Balogh et al. (2001) derive an 
even lower value of about 5%, based on the A"-band lumi- 
nosity function (Cole et al., 2001), which should in principle 
provide a more robust estimate due to the relatively weak 
dependence of the K-band mass-to-light ratio on star for- 
mation history. However, some uncertainty arises from the 
adopted IMF; as Cole et al. (2001) discuss, an estimate of 
the mass fraction in stars higher by about a factor of 2 is 
obtained if a Salpeter IMF (1955) instead of a Kennicutt 
(1983) IMF is assumed. Note that there are also tight limits 
on the amount of gas in neutral and molecular form, imply- 
ing that they together can only account for at most 10% of 
the mass bound in stars. Hence, most of the condensed gas 
should in fact be bound in stars, and there is not much room 
to "hide" gas that has cooled. 

The discrepancy between the high efficiency of cooling 



and the low abundance of stars suggests that any feedback 
model lacking a mechanism for returning baryons from the 
condensed phase into the extended halo or IGM will be un- 
able to account for the observations. For example, in the no- 
wind runs shown in Fig. 13, more than 18.2% of the baryons 
are locked up in stars, which is already above current ob- 
servational constraints, even though this number represents 
an underestimate due to the limited resolution of these sim- 
ulations, and it does not yet include the dense, neutral gas 
remaining in the simulation. 

However, winds are able to significantly affect the cos- 
mic star formation history. In Figure 14, we show simulation 
results for the comoving star formation density as a func- 
tion of redshift for different wind strengths. The dotted line 
shows a 2 x 50 3 run using our fiducial set of parameters and 
the quiescent model of star formation without winds. The 
solid line is the same model, but this time including winds 
with a strength specified by r\ = 2 and x = 0.25, i.e. star- 
forming regions blow a wind with speed 242 km s -1 at a 
mass-loss rate equal to twice the star formation rate. This 
reduces the integrated star formation rate such that 13.5% 
of the baryons become bound in stars, providing a better 
match to direct observational constraints of the cosmic lu- 
minosity density, although perhaps still being uncomfortably 
high. 

Faster winds can escape from bigger haloes and provide 
larger heating to the surrounding IGM. This can make them 
more effective in reducing subsequent star formation, even 
if the mass-loss rate is smaller. For example, the dashed line 
in Fig. 14 shows a rather extreme wind model, with r\ = 0.5 
and x = 1- The wind has thus an initial speed of 968 kms -1 , 
but a mass-loss rate four times smaller than in our fiducial 
model. This energetic wind leads to substantial heating of 
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the IGM, and reduces the star formation rate considerably, 
such that this model turns only 8% of its baryons into stars. 

Finally, we also show the result for a 2 x 100 3 simulation 
of our fiducial model. Here, the star formation rate is sub- 
stantially larger at high redshift. This is expected because 
the improved resolution now allows much more of the star 
formation to be seen in low-mass objects that begin to form 
abundantly at high redshift. It is thus clear that much better 
resolution than was used in the 2 x 50 3 runs is required to 
obtain converged results in the high-redshift regime. How- 
ever, this is not necessarily the case for the integrated star 
formation rate. Since the elapsed time is rather small at 
high redshift, the amount of stellar material formed then is 
small compared to the total, despite the large star formation 
rates. Therefore, the 2 x 100 3 run produces only moderately 
more stars by z = 0; the fraction of baryons is 14.8%, up 
from 13.5% for the corresponding 2 x 50 3 run. Note that the 
slight upturn seen in the 2 x 50 3 and 2 x 100 3 runs with slow 
winds at redshifts below ~ 0.5 is due to the small box size 
of these simulations, which do not allow a proper sampling 
of the halo mass function at low redshift, where in fact the 
fundamental modes of the boxes already become non-linear. 
In particular, the star formation rate begins to be domi- 
nated by the few most massive halos in the box, and they 
are able to reaccrete at low redshift some of the gas ejected 
by winds out of their progenitor halos. If larger cosmological 
volumes are simulated, such biases in the estimated mean of 
the cosmic star formation density can be avoided (Springel 
& Hernquist, 2002b). 



6.3 Metal enrichment of the IGM 

As we discussed in Section 4, galactic winds may be of crucial 
importance for the transport of metals into the IGM, where 
they have been observed in absorption line systems down to 
very low column densities. In Figure 15, we show the mean 
metallicity of the g cLS clS cl function of baryonic overdensity 
in our cosmological simulations. We give results at redshifts 
z — 4, 2 = 2.33, and z = 0, for runs carried out with our 
fiducial parameters, both with and without winds. 

It is evident that in the simulation without winds, met- 
als are nearly confined to gas with overdensities S > 10. 
At redshift z = 0, such gas reaches a metallicity of about 
1O~ 4 Z , where Z Q is the solar metallicity. However, at red- 
shifts higher than z > 2.3 even gas at overdensities 100 
is well below 10~ 4 Z Q , and thus falls significantly short of 
the metallicities that are observed in the Ly-a forest at 
these redshifts, as indicated by the shaded region (following 
Aguirre et al., 2001a). Recall that in our self-regulated model 
for star formation, stars begin to form only at high physical 
densities, corresponding to baryonic overdensities of about 
~ 10 6 at z = 0. If winds are not included, only dynamical 
stripping may bring metal-enriched gas from the ISM to the 
lower density environments. The results of Fig. 15 suggest 
that these processes are not efficient enough to explain the 
enrichment of the IGM. 

However, including winds changes the enrichment pat- 
tern significantly, as expected. Our fiducial model enriches 
the IGM to an interesting metallicity of Z ~ 10~ 2,5 Zq for 
the relevant densities. It is thus clear that a wind model like 
the one discussed here can in principle account for the mean 
metallicity of the IGM. 
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Figure 15. Mean metallicity of the g function of baryonic 

overdensity. The three panels show results for redshifts 2 = 4, 
2 = 2.3, and 2 = 0, both for simulations with and without winds. 
We follow Aguirre et al. (2001a) and use a shaded region to ap- 
proximativcly indicate the range of metallicities observed in Ly-a 
absorption line studies at z ~ 3. 



Note that the distribution of metals in the gas is highly 
inhomogeneous. This is seen in Figure 16, where we show 
projected metallicity maps at redshifts z = 2.3 and 2 = 0. 
The detailed metallicity distribution together with the in- 
homogeneous heating pattern due to the winds potentially 
yield signatures in the Ly-a forest that may be very con- 
straining for the enrichment model discussed here. We plan 
to address this question further in future work. Note that at 
low velocities, winds will primarily be able to escape from 
small haloes, pushing the epoch of enrichment of the IGM 
to high redshift, when these haloes form abundantly, and 
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Figure 16. Projected mean metallicity of the gas in a 2 X 50 3 simulation that includes galactic winds. The map is 11.3 h~ x Mpc on a side 
and shows the full simulation box in projection at rcdshift z = 2.3 (left) and z = (right). The mean metallicity of the gas is indicated 
by a logarithmic colour-scale. 



when the winds can propagate relatively far because of the 
small scale-size of the universe at the epoch of ejection. 



7 DISCUSSION 

We have presented a new model for the treatment of star for- 
mation and feedback in cosmological simulations of galaxy 
formation. Our approach leads to the establishment of a 
tight self-regulation cycle for star formation which is based 
on a rough, yet physically motivated model of the ISM. A 
crucial ingredient of the model is the assumed multi-phase 
structure of the ISM. At high density, much of the interstel- 
lar material becomes bound in cold clouds, becoming sites 
for star formation. The density of the ambient phase is low- 
ered accordingly, which also reduces its radiative losses and 
allows it to be heated by supernovae to a temperature high 
enough to provide some pressure support for the ISM. The 
stabilising effect of this pressure counteracts star formation 
and, together with the cloud formation/evaporation cycle, 
leads to self-regulation. Due to the nature of the approxi- 
mations made, it is clear that our model is still phenomeno- 
logical to a large degree. In principle, however, the physical 
approximations we used can be refined in the future, thereby 
improving the faithfulness of the model. 

The star formation timescale in the quiescent mode of 
star formation can be directly determined from observations 
of local disk galaxies. If this normalisation is adopted, we 
find that cosmological simulations nevertheless lead to an 
overproduction of stars. We have invoked galactic winds as 
a heuristic extension of our model to remedy this problem. 
While there is mounting observational evidence for the ex- 
istence of such winds, they are introduced in our model in 



a phenomenological manner. However, winds clearly have a 
range of highly interesting consequences. They reduce the 
global efficiency of star formation, and this suppression is 
particularly strong in low-mass haloes, thereby helping to 
explain the observed "low" values of the luminosity density. 
Winds are also capable of accounting for the enrichment of 
the low-density IGM with metals, and they may be crucial 
to understanding the metal distribution in the intra-cluster 
gas of clusters of galaxies. 

Another important property of our technique is that it 
is numerically well posed in the sense that the parameters 
of the model can be determined directly based on physi- 
cal arguments or observational input, and need not to be 
changed when the mass resolution is varied. The use of the 
sub-resolution technique introduces a well-specified contin- 
uum limit in the hydrodynamic equations, despite the inclu- 
sion of cooling. Unlike in standard single-phase simulations, 
an unphysical collapse of gas to a scale set by the resolution 
limit of the simulation is thus prevented. Combined with im- 
proved formulations of SPH (Springel & Hernquist, 2002a) 
this makes it much easier to reach numerical convergence 
even under conditions of only moderate resolution. 

It will thus be interesting to study the predictions of our 
model in more detail. As first steps, we have already carried 
out a detailed analysis of the predictions of the model for 
the cosmic star formation history in the ACDM cosmology 
(Springel & Hernquist, 2002b), and we analysed the im- 
pact of cooling, star formation, and winds on secondary 
anisotropies of the cosmic microwave background (White 
et al., 2002). In these studies, a new comprehensive set of 
simulations that included strong galactic winds was used. 
We plan to extend our analysis of these simulations in fu- 
ture work. 
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